← LAB ◊ ↑ TOC ↑

Name, Vorname:
Hofmann Fabian
Richter Daniel

LAB 3: Filterentwurf mit pyfda

In diesem Versuch

Inhalt

0. Allgemeine Hinweise

1. Überblick über pyfda

2. Filterentwurf

3. Filterexport und -simulation

0. Allgemeine Hinweise

0.1 Jupyter Notebook

• \<SHIFT>-\<RETURN> führt eine Codezelle aus und rendert eine Textzelle.
• In Markdown sind Leerzeilen wichtig zum Trennen von Abschnitten!
• Sie können LaTeX-Code zwischen \$ ... \\$ einschließen.
• Kontexthilfe zu Funktionen etc. bekommen Sie über \<SHIFT>-\<TAB>

Abbildungen in diesem Notebook wurden konvertiert mit https://www.base64-image.de/ und in den HTML-Code eingebettet.

0.2 Installation von pyfda

Für diesen Versuch benötigen Sie die Software pyfda (Python Filter Design and Analysis), die Sie per pip install pyfda installieren können (wenn Sie eine Python-Installation auf Ihrem Rechner haben), alternativ gibt es unter https://github.com/chipmuenk/pyfda Binaries für Windows und für Ubuntu 20.10 (funktioniert vermutlich auch mit anderen halbwegs aktuellen Distros).

0.3 Abgabe

Nach dem Praktikumsversuch exportieren Sie das Notebook mit Textantworten, Codezellen und Plots als HTML (File -> Export Notebook As ... -> Export Notebook to HTML) und reichen es in Moodle ein.</div>

---

import os, sys
module_path = os.path.abspath(os.path.join('..')) # append directory one level up to import path
if module_path not in sys.path:  # ... if it hasn't been appended already
    sys.path.append(module_path)
import dsp_fpga_lib as dsp
dsp.versions() # print versions 
    
%matplotlib inline
import matplotlib.pyplot as plt
size = {"figsize":(12,5)} # Plotgröße in Inch

import numpy as np
import scipy.signal as sig
import wave

from IPython.display import Audio, display
#-----------------------------------------------------------------------------
def wav2np(filename):
    """ Read the wav-file and convert it to a one or two-dimensional numpy array,
        depending on the number of channels. 
        Properties of the WAV-file are stored as function attributes (evil)
    """
    wf = wave.open(filename,'rb')
    wav2np.N_CH = wf.getnchannels() # number of channels
    wav2np.W    = wf.getsampwidth() # wordlength per sample in bytes
    wav2np.N_FR = wf.getnframes()   # number of frames
    wav2np.f_S  = wf.getframerate() # sample (frame) rate 
    
    print("{0} channels with {1} frames of {2} bytes and f_S = {3} Hz.".format(wav2np.N_CH, wav2np.N_FR, wav2np.W, wav2np.f_S))
    
    if wav2np.W == 2:
        samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int16) # read wav data as 16 bit integers, R and L samples are interleaved
    elif wav2np.W == 1:
        samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int8) # read wav data as 8 bit integers, R and L samples are interleaved
    else:
         raise TypeError("Unknown data format: {0} bytes".format(wav2np.W))   

    samples = np.array([samples_in[idx::wav2np.N_CH] for idx in range(wav2np.N_CH)], dtype=np.int32) # deinterleave channels to numpy array N_CHAN x N_FRAMES
    return samples

Python version: 3.10.8
Numpy: 1.23.5
Scipy: 1.9.3
Matplotlib: 3.6.2 module://matplotlib_inline.backend_inline

1. Überblick über pyfda

Auf der linken Seite wählen und entwerfen Sie Filter unter den Tabs "Specs" (Filterentwurf), "b,a" (Eingabe von Filterkoeffizienten), "P/Z" (Eingabe von Polen und Nullstellen). Auf der rechten Seite schauen Sie sich Betrags- und Phasengang, Gruppenlaufzeit, P/N Diagramm und die transiente Antwort auf verschiedene Stimuli an.

1.1. Frequenzen und Amplituden

Frequenzen können in unterschiedlicher Form eingeben und angezeigt werden:

• Normalisiert auf die Samplingfrequenz $f_S$: $F = f/f_S$
• Normalisiert auf die Nyquistfrequenz $f_{Ny} = f_S/2$: $F = 2f/f_S$
• In absoluten Frequenzen, hier können Sie die Abtastfrequenz $f_S$ eingeben. In der Einstellung "Lock" (verriegeltes Schloss) bleiben die absoluten Frequenzen gleich wenn die Abtastrate geändert wird. Bei "Unlocked" werden die absoluten Frequenzen so angepasst, dass die normalisierten Frequenzen gleich bleiben und damit auch Filterentwürfe, Signalverläufe etc.

Amplituden werden in V, W oder dB eingegeben und angezeigt.

1.2 Filterspezifikationen

Filterentwurfsalgorithmen benötigen unterschiedliche Eingabeparameter. Wenn "Order Minimum" angeklickt ist, versucht ein Algorithmus die minimale Ordnung zu den gegebenen Filterspezifikationen zu berechnen, halbwegs genau funktioniert das aber nur bei Equiripple (FIR) und den meisten IIR-Filtern. Bei pyfda können Sie nach dem Entwurf mit minimaler Ordung das Häkchen entfernen und dann selber Entwurfsparameter bearbeiten.

Nicht benötigte Parameter sind entweder deaktiviert oder ausgeblendet. Manche Parameter sind ausgegraut, aber bearbeitbar: Diese Parameter beeinflussen nicht den Filterentwurf, sie dienen aber z.B. dazu um das Toleranzschema anzuzeigen.

2. Filterentwurf

Ein Audiotrack von einer CD ($f_S = 44,1 \,\text{kHz}$) soll tiefpassgefiltert werden. Die Filterspezifikationen sind:

• Eckfrequenz des Passbands: $f_{PB} = 3,4\, \text{kHz}$
• Eckfrequenz des Sperrbands: $f_{SB} = 4\, \text{kHz}$
• Max. Ripple im Passband: $A_{PB} = 1\, \text{dB}$
• Min. Dämpfung im Sperrband: $A_{SB} = 60\, \text{dB}$

Es sollen die folgenden Filtertypen getestet werden:

• Fourier-Approximation / Least-Square (Windowed FIR mit Boxcar-Fenster)
• Windowed FIR mit Kaiser-Fenster
• Equiripple
• Butterworth
• Chebychev 1
• Chebychev 2
• Elliptisch (Cauer)

2.1 Filtervergleich

SIMULATION:

Ermitteln Sie die folgenden Filtereigenschaften und fassen Sie sie in einer Tabelle zusammen mit den Spalten:

Filtertyp	IIR/FIR?	Ordnung N	Gruppenlaufzeit τg	Mult. L pro Sample	Max. Samplerate [GHz]	Vorteil	Nachteile
Fourier-Approximation / Least-Square(Windowed FIR mit Boxcar-Fenster)	FIR	164	1.84 ps	164	1.83	geringe Nichtlinearität	hohe Rippel / zu geringe Steilheit
Windowed FIR mit Kaiser-Fenster	FIR	268	30.2 ns	268	1.12	keine Rippel im Sperrbereich / hohe Däpfung im Sperrbereich	mittlerer Rippel im Sperrbereich
Equiripple	FIR	163	18.4 ps	163	1.8	geringe Nichtlinearität	hoher Rippel
Butterworth	IIR	PyFDA -> max. 20 f_SB = 4.7kHz	1.46 ms	23	13.0	kein Rippel / geringe Nichttlinearität	geringe Steilheit
Chebychev 1	IIR	14	4.75 ms	17	17.6	mittlere Steilheit	hoher Rippel im Durchlassband / hohe Nichttlinearität
Chebychev 2	IIR	14	1.48 ms	17	17.6	mittlere Steilheit	hoher Rippel im Sperrbereich / hohe Nichttlinearität
Elliptisch (Cauer)	IIR	7	2.31 ms	10	30	geringe Ordung	hoher Rippel

• Welche der Filter sind FIR, welche IIR?

• Vergleichen Sie die benötigte Ordnung N. Für nicht alle der Filter gibt es gut funktionierende Designalgorithmen, Sie müssen dann von Hand die Filterordnung und / oder Eckfrequenzen so lange anpassen, bis alle Spezifikationen erfüllt sind. Im Tab "Info" sehen Sie, ob Spezifikationen verletzt werden.

• Welche Ordnung benötigen Sie, wenn Sie die Abtastrate vor der Filterung auf die Hälfte reduzieren? (Tipp: Aktivieren Sie das "Schloss" neben der Samplingfrequenz, bevor Sie f_S ändern.) Es genügt, wenn Sie das für ein Filter ausprobieren. Was müssen Sie beachten wenn Sie die Samplerate reduzieren wollen?

• Welche Gruppenlaufzeit τ_g haben die Filter bei $f = 2 \,\text{kHz}$? Schauen Sie sich im y[n] an, wie die Signalform eines 500 Hz Sägezahnsignals durch das Filter verformt wird.

• Wieviele Multiplikationen L pro Sample werden benötigt? Berücksichtigen Sie dabei die Filterstruktur (Optimierung durch Linearphasigkeit möglich?), IIR Filter sollen in einer kaskadierten SOS-Struktur realisiert werden - siehe folgende Abbildungen.

  • FIR: L = N
  • IIR L = N + 3
    
    

• Welche maximale Samplerate kann man auf einem DSP mit max. $300 \cdot 10^6$ Multiplikationen / s verarbeiten? Zusätzliche Informationen: How fast are DSPs?

  • f_s_max = f_DSP / L

• Fügen Sie den Plot des Betragsgangs mit Phasengang (Checkbox "Phase") ein.
  
  
  
  
  
  Windowed FIR mit Boxcar-Fenster:
  
  
  
  
  
  Windowed FIR mit Kaiser-Fenster:
  
  
  
  
  
  Equiripple:
  
  
  
  
  
  Butterworth:
  
  
  
  
  
  Chebychev 1:
  
  
  
  
  
  Chebychev 2:
  
  
  
  
  
  Elliptisch:
  
  

Direktform	Lin. phasige Form	Kaskadierte Second Order Sections

3. Filterexport und -simulation

Im folgenden sollen Sie Audiofiles mit dem Equiripple und dem elliptischen Filter filtern. Dazu können Sie aus pyfda Filter in verschiedenen Formaten exportieren:

• Im Tab Specs speichern Sie Filter komplett ab, um sie später wieder in pyfda zu laden
• Im Tab b,a exportieren Sie die Filterkoeffizienten, u.a. im CSV (Comma Separated Values) Format. Sie können mit dem "Zahnrad" einstellen, dass die Koeffizienten direkt zum Clipboard (Zwischenablage) exportiert werden. Von dort können Sie sie einfach als Python Listen oder Arrays (a = [...]) in die nächste Codezelle einfügen.
• Im Tab P / Z exportieren Sie die Pole und Nullstellen, auch hier können Sie über die Zwischenablage gehen.

Sie können Daten im Zeitbereich filtern mit

y = sig.lfilter(b,a,x)

b und a sind die Koeffizientenvektoren (nicht-rekursiver / rekursiver Teil), x ist das Eingangssignal. y hat die gleichen Dimensionen wie x (es können auch Stereofiles verarbeitet werden). Was müssen Sie bei einem FIR-Filter für a angeben?

* a = [1, 0, 0, ....] = [1]

Wenn die Filterbeschreibung als Second-Order Sections vorliegt, kann auch der folgende Befehl benutzt werden (allerdings unterstützt pyfda noch keinen Export im sos - Format):

y = sig.sosfilt(sos,x)

Anhören können Sie sich die gefilterten Audiofiles wieder mit display(Audio(data=y, rate=wav2np.f_S)).

Sie können auch eigene WAV-Files in Ihr medien Verzeichnis hochladen und testen, bei großen Audiofiles sollten Sie zunächst nur ein paar Sekunden testhören (mit x[:N_samples] bzw. x[:, :N_samples]) bei Stereo-Samples. Was bewirkt die seltsame Syntax bei Stereosignalen?

Hören Sie einen Unterschied zwischen dem Orginal und den beiden unterschiedlichen Filtern? Woran liegt das?

• Orginal vs. Filter -> gefiltert Dupfer (Frequenzen über 4kHz um >60dB gedäpft)

x_w = wav2np("../medien/87778__marcgascon7__vocals.wav")

#Equiripple
b_FIR = [0.0006123033469102477,0.0005231312265310316,0.000561083771645731,0.00041518098768839355,2.827711686114232e-05,-0.0006169795177762187,-0.0014804911786779923,-0.002458026960357289,-0.0033955247678761798,-0.004109019417673139,-0.004428667554745472,-0.004237079147342417,-0.0035082858849691473,-0.002325963428206755,-0.0008793225682988921,0.0005669966191960529,0.001727314146414968,0.002358436650440632,0.0023225886047265588,0.0016281694762668173,0.00044002268486894485,-0.0009495414647988306,-0.00218391145568962,-0.002926600347304991,-0.002948523144002793,-0.0021936643092102604,-0.0008066063399116822,0.0008911226757635836,0.002471991422020751,0.0035073292743326347,0.0036806624276502827,0.0028800941086848594,0.0012444616340492698,-0.0008540601244899048,-0.0028900585715064993,-0.0043130326782202496,-0.004692235633235503,-0.00384162591921165,-0.0018896284603115599,0.0007297360484630355,0.0033677431056925913,0.005317244868695226,0.005996942088874641,0.005117766406001423,0.0027847237841484847,-0.0004979172934868248,-0.0039276848931389715,-0.0065934641947688385,-0.0077099517214559606,-0.006839091317099379,-0.004037685130741728,0.00011364684868242633,0.004618634275870731,0.008293677454136077,0.010066529782336353,0.009273031150922878,0.005873052144054021,0.0005221554755543486,-0.005533741821921171,-0.010729128825649713,-0.013560660832545542,-0.012991437230981336,-0.008775651715744494,-0.00161273198833575,0.0069258868738092455,0.014701407248243398,0.01949170936897492,0.019552878251152542,0.014129162521417389,0.003783120765910111,-0.00955025960828738,-0.02282413158635936,-0.03242421406707713,-0.03489768645270763,-0.027719185939259615,-0.009927599299774674,0.017513404923713328,0.051725839878917296,0.08831767123202569,0.12210814306160755,0.14804243901826264,0.1621149051862032,0.1621149051862032,0.14804243901826264,0.12210814306160755,0.08831767123202569,0.051725839878917296,0.017513404923713328,-0.009927599299774674,-0.027719185939259615,-0.03489768645270763,-0.03242421406707713,-0.02282413158635936,-0.00955025960828738,0.003783120765910111,0.014129162521417389,0.019552878251152542,0.01949170936897492,0.014701407248243398,0.0069258868738092455,-0.00161273198833575,-0.008775651715744494,-0.012991437230981336,-0.013560660832545542,-0.010729128825649713,-0.005533741821921171,0.0005221554755543486,0.005873052144054021,0.009273031150922878,0.010066529782336353,0.008293677454136077,0.004618634275870731,0.00011364684868242633,-0.004037685130741728,-0.006839091317099379,-0.0077099517214559606,-0.0065934641947688385,-0.0039276848931389715,-0.0004979172934868248,0.0027847237841484847,0.005117766406001423,0.005996942088874641,0.005317244868695226,0.0033677431056925913,0.0007297360484630355,-0.0018896284603115599,-0.00384162591921165,-0.004692235633235503,-0.0043130326782202496,-0.0028900585715064993,-0.0008540601244899048,0.0012444616340492698,0.0028800941086848594,0.0036806624276502827,0.0035073292743326347,0.002471991422020751,0.0008911226757635836,-0.0008066063399116822,-0.0021936643092102604,-0.002948523144002793,-0.002926600347304991,-0.00218391145568962,-0.0009495414647988306,0.00044002268486894485,0.0016281694762668173,0.0023225886047265588,0.002358436650440632,0.001727314146414968,0.0005669966191960529,-0.0008793225682988921,-0.002325963428206755,-0.0035082858849691473,-0.004237079147342417,-0.004428667554745472,-0.004109019417673139,-0.0033955247678761798,-0.002458026960357289,-0.0014804911786779923,-0.0006169795177762187,2.827711686114232e-05,0.00041518098768839355,0.000561083771645731,0.0005231312265310316,0.0006123033469102477]
a_FIR = [1.]


#Elliptisch
b_IIR = [0.001710952004172366,-0.005755838859144107,0.008363025285148405,-0.004112746893907044,-0.004112746893907043,0.008363025285148403,-0.005755838859144107,0.001710952004172366]
a_IIR = [1.0,-6.10402966946488,16.40268768530511,-25.100822235115835,23.595177173520497,-13.615049072423911,4.4641836431418325,-0.6417367418902715]


#Signal Filtern
y_FIR = sig.lfilter(b_FIR, a_FIR, x_w)
y_IIR = sig.lfilter(b_IIR, a_IIR, x_w)

print("Orginal")
display(Audio(data=x_w, rate=wav2np.f_S))


print("-----------------------------------------------------------")
print("FIR: Equiripple")
print(y_FIR.dtype)
display(Audio(data=y_FIR, rate=wav2np.f_S))


print("-----------------------------------------------------------")
print("IIR: Elliptisch")
print(y_IIR.dtype)
display(Audio(data=y_IIR, rate=wav2np.f_S))

2 channels with 349952 frames of 2 bytes and f_S = 44100 Hz.
Orginal

-----------------------------------------------------------
FIR: Equiripple
float64

-----------------------------------------------------------
IIR: Elliptisch
float64

#fig, ax = plt.subplots(**size)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(35,8))
l = 4000
start = 20000


ax[0].set_title("FIR-Filter: Equiripple")
ax[0].plot(np.arange(l), x_w[0][start:start+l], label="Orginal")
ax[0].plot(np.arange(l), y_FIR[0][start:start+l], label="y_FIR")
ax[0].legend()

ax[1].set_title("IIR-Filter: Elliptisch")
ax[1].plot(np.arange(l), x_w[0][start:start+l], label="Orginal")
ax[1].plot(np.arange(l), y_IIR[0][start:start+l], label="y_IIR")
ax[1].legend()

<matplotlib.legend.Legend at 0x24d37d0a260>

---

Copyright

(c) 2016 - 2021 Prof. Dr. Christian Münker

This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. The latest version can be found at https://github.com/chipmuenk/dsp.

This notebook is provided as Open Educational Resource. Feel free to use it for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Christian Münker, Digital Signal Processing - Vorlesungsunterlagen mit Simulationsbeispielen, 2020.